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Abstract 

In the case of ground-based telescopes equipped with adaptive optics systems, the point spread function (PSF) is only poorly 
known or completely unknown. Moreover, an accurate modeling of the PSF is in general not available. Therefore in several 
imaging situations the so-called blind deconvolution methods, aiming at estimating both the scientific target and the PSF from the 
detected image, can be useful. A blind deconvolution problem is severely ill-posed and, in order to reduce the extremely large 
number of possible solutions, it is necessary to introduce sensible constraints on both the scientific target and the PSF. 

In a previous paper we proposed a sound mathematical approach based on a suitable inexact alternating minimization strategy 
for minimizing the generalized Kullback-Leibler divergence, assuring global convergence. In the framework of this method we 
showed that an important constraint on the PSF is the upper bound which can be derived from the knowledge of its Strehl ratio. 
The efficacy of the approach was demonstrated by means of numerical simulations. 

In this paper, besides improving the previous approach by the use of a further constraint on the unknown scientific target, we 
extend it to the case of multiple images of the same target obtained with different PSFs. The main application we have in mind 
is to Fizeau interferometry. As it is known this is a special feature of the Large Binocular Telescope (LBT). Of the two expected 
interferometers for LBT, one, LINC-NIRVANA, is forthcoming while the other, LBTI, is already operating and has provided the 
first Fizeau images, demonstrating the possibility of reaching the resolution of a 22.8 m telescope. Therefore the extension of our 
blind method to this imaging modality seems to be timely. 

The method is applied to realistic simulations of imaging both by single mirrors and Fizeau interferometers. Successes and 
failures of the method in the imaging of stellar fields are demonstrated in simple cases. These preliminary results look promising 
at least in specific situations. The IDL code of the proposed method is available on request and will be included in the forthcoming 
version of the Software Package AIRY (v.6.1). 
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1. Introduction 


In a previous paper ( Prato et akl 2013 ) we proposed a blind 
deconvolution method for ground-based telescopes equipped 
with an adaptive optics (AO) system. Assuming that the image 
and the corresponding background are known, then the features 
of the method are the following: 


• formulation of the problem as a constrained minimization 
of the data fidelity function in the case of Poisson noise 
(photon counting), namely a generalized Kullback-Leibler 
divergence depending on the unknown astronomical target 
(in the following called the object) and on the unknown 
point spread function (PSF); 


• non-negativity of the object; 
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non-negativity of the PSF and normalization to unit vol¬ 
ume; 


upper bound on the PSF values derived from the 
knowledge of the Strehl ra tio (SR), as suggested by 


Desidera and CarbilleJ ( 200^ : 


use of th e inexact alternat ing minimization method pro¬ 
posed by Bonettinil ( 201 ih and ba sed on the scal e d gra¬ 


dient projection (SGP) method ( Bonettini et al. . 20091 : 


Prato et alll2012l) . 


The method is iterative and at each outer iteration the object and 
the PSF are updated by means of given (but arbitrary) numbers 
of inner iterations of SGP. We remark that when SGP is applied 
to the object only projection on the non-negative orthant is re¬ 
quired while, in the case of the PSF, the projection is performed 
on the convex and closed set defined by the box and equality 
constraints. 

Our proposed method is similar to bli nd methods based 
on the Rich ardson-Lucy (RL) algorithm ( Richards'onl 1972 : 
Lucy , 19741) or its accelerated versions, suc h as those p r opose d 
in Holmes ( 1992 ): Tsumurava et al. ( 19941) : Fish et al. (1995); 
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Biggs and AndrewsI ( 19981) for the case of single image and in 
Desidera et alJ ( 2006 ) for the case of multiple images, applica¬ 


ble to Fizeau interferometry. But the advant age of usi ng SGP 


ant age ot usi ng , 
in Bonettinil (1^ 


in place of RL is double; first, as proved in [Bonettinil (1201 11) . 
global convergence of the iteration holds true, i.e. any limit 
point of the sequence is a stationary point of the constrained KL 
divergence; secondly it is possible to introduce box and equal¬ 
ity constraints on the different blocks of variables (object and 
PSF). 

Therefore the novelty of the method is that it is based on a 
sound mathematical approach and allows, in an easy way, the 
introduction of important constraints on the PSF such as the 
Strehl constraint, preventing the appearance of trivial solutions 
such as a delta function. 

The method, in its present form, does not consider the use 
of regularization for the object or the PSF, in addition to that 
provided by the constraints mentioned above. We remark that 
it is quite easy to introduce as an additional constraint on the 
object the value of its flux (namely, its (\ norm), as derived from 
the detected images and the knowledge of the background. This 
constraint, enforcing the sparsity of the object, is considered in 
this paper. If we select a large number of inner iterations on the 
object variables, the method is suitabl e for the reconstru ction 
of star systems, as already remarked in lPrato et al.l (l2013h : this 
result is confirmed in this paper by means of simulations more 
realistic than those used in that paper. 

However the main contribution of this paper is the extension 
of the method to the case of Fizeau interferometry. As it is 
known this is a specific feature of the Large Binocular Tele¬ 
scope which consists of two 8.4 m mirrors situated on a com¬ 
mon mount with a center to center distance of 14.4 m. Indeed, 
this structure is suitable for Fizeau interferometry which should 
provide images with the resolution of a 22.8 m telescope in the 
direction of the baseline joining the center of the two mirrors 
and that of a 8.4 m telescope in the orthogonal direction. Dif¬ 
ferent images of the same target corresponding to different ori¬ 
entations of the baseline can be combined by suitable decon¬ 
volution methods to provide a unique r econstructed image w ith 


the highest resolution in all directions (iBertero et al.Ll201 ll) . 


Two interferomete rs are planned for LBT: the forthcom¬ 
ing LINC-NIRVANA ( Herbst et al. . 2003 ). in advanced realiza¬ 
tion stage by a German-Italian consortiu m leaded by MPIA , 


Heidelberg, and th e NASA funded LBTI (IWilson et al.L 12008 


Bailey et~^ 2014 ) already operating on Mount Graham. In¬ 
deed, images of the Jupiter moon lo were obtained with 
LBTI/LMIRcam during UT 2013 Dec ember 24, showing tha t 


the resolution of a 22.8 m is reachable (iLeisenring et al.Ll2014h 


and thus proving that LBT is the first in a class of extremely 
large telescopes (ELT). 

All the methods developed for Fizeau interferometry are also 
applicable to other situations where multiple images of the 
same target, corresponding to different PSFs, are available such 
as the co-adding problem in Astronomy ( Lucv and Hood. 1992 ) 
or the multiple image method use d in STEP microscop y for im¬ 
proving the signal-to-noise ratio (ICastello et al.L 120141) . 


Eor this reason we present in Sect. 2 our blind method with¬ 
out a specific reference to Eizeau interferometry but just as a 


method for multiple image deconvolution in the case of Pois¬ 
son noise, including, as it is obvious, the case of a single image 
as a particular case. In the same section we discuss the intrinsic 
limitations of our constrained blind deconvolution, a discussion 
which is possible in our sound mathematical framework. 

In the simulations intended to validate the method we focus 
on LBT which is equipped with a very innovat ive AO system. 


the so -called Eirst Light AO (ELAO) system (lEsnosito et al 


20 id) , providing SR values up to 0.9 in K-band. Therefore for 


single image simulations we use models of the PSE of such a 
system. On the other hand, for multiple image deconvolution 
we consider images generated by means of PSEs computed for 
the interferometer LINC-NIRVANA (LN) in K band. Since the 
camera of LINC-NIRVANA has a pixel size of ~ 5 mas while 
the LMIRcam of LBTI has a pixel size of 10.7 mas, the shape 
of the PSEs of LN in K band (2.2/im) is similar to the shape of 
the PSEs of LBTI in M band (4.8/im). Obviously the properties 
of the images may be very different. Details on image modeling 
and simulation are given in Sect. 3. 

Einally in Sect. 4 we discuss our numerical results in the case 
of binary systems and “open cluster” models. Conclusions are 
sketched in Sect. 5. 


2. Method 


We assume that p different images of the scientific object, 
with p different PSPs are available. The case of a single aper¬ 
ture telescope obviously corresponds to p = 1 if only one im¬ 
age has been acquired; if different observations have been per¬ 
formed at different times, hence with PSPs corresponding to 
different AO corrections, then the approach can be used for the 
co-adding of these images. 

Let / be the unknown astronomical object and let Kj be the 
unknown PSPs (each one normalized to unit volume) corre¬ 
sponding to the detected images gj for j - l,...,p (we as¬ 
sume for simplicity a space-invariant model), then we define 
as Ajf - Kj*f the corresponding imaging matrices. Moreover 
we denote as b j the expected value of the background emission 
in image gj and we assume that it is known, so that the expected 
value of gj is given by Ajf + b j. 

Since it is quite natural to assume that the p images are statis¬ 
tically independent, the likelihood of the problem is the product 
of the likelihoods of the different images. We assume that they 
are perturbed by Poisson noise. Then, by taking the negative 
logarithm of the likelihood we obtain the following data-fidelity 
function which is the sum of p Kullback-Leibler (KL) general¬ 
ized d ivergences, also known as Csiszar I-divergences (iCsiszar , 
19911) . one for each image, i.e. 


Jo(f,Ki,...,Kp-g,b)^ 

’’ gj(m) 


(A/)(m)H-ii/m)-g/m)} , 


( 1 ) 


where S is the set of the values of the multi-index m character¬ 
izing the pixels of the image array, and ig,b) = {{gj,bj)yj^^. 
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The problem of image deconvolution (without regularization) 
consists in the minimization of this function with respect to 
/ for given PSFs, images and backgrounds. The minimiza¬ 
tion can be obtained by means o f RL in the single image 
case, by means of OS - EM m ethod dHudson and Larkin . 1994 : 
Bertero and Boccaccil 2000) in the multiple image case, or 
by rn eans of SGP method ( Bonettini et all 20091: Prato et al 
2012 ) in both cases. As shown in Prato et al. ( 20121) SGP i 


more efficient than OS-EM if the number p of images is not too 
large. 

2.7. Blind deconvolution: problem formulation 

If images and backgrounds are given, the ML approach to 
blind deconvolution can be formulated as the minimization of 
the function in Eq. dB with respect to p -i-1 blocks of unknown 
variables, namely the object / and the p PSEs Kj, j = 1, p. 
As it is known, this function is convex with respect to each 
block of variables for fixed values of the ot hers, but is not con - 
vex with respect to the full set of variables (IPrato et al.Ll2013l) . 
Therefore blind deconvolution is a difficult problem of noncon- 
vex optimization. Moreover, this problem is highly ill-posed 
and allows uninteresting solutions. Eor instance, a global min¬ 
imum can be achieved by choosing / - g — b and Kj - 6, 
j - 1,..., p, where 6 is the Dirac delta array. This trivial mini- 
mizer can be avoided by the introduction of suitable regulariza¬ 
tion terms and constraints. 

Since we mainly consider the case of stellar fields or, in other 
words, of sparse objects, by taking into account the sparsity 
property of the mini mizers of the KL dive rgence in the case of 
image deconvolution ( Bertero et all 2009 ). we do not introduce 
an object-dependent regularization term in the objective func¬ 
tion. However, besides non-negativity of the object we also in¬ 
troduce a constraint on its flux; more precisely we require that 
the object flux coincides with the average flux of the p detected 
images (after background subtraction), which is given by 




( 2 ) 


j=\ mG5 


We remark that this consttaint is further enforcing sparsity; in 
the case of deconvolution and zero value of the backgrounds, 
it is automatically satisfied by the minimizers of the KL diver¬ 
gence. 


As conc erns the PSEs, as shown in iDesidera and Carbillel 
( 2009h and Prato et al. ( 2013 ). an important constraint is the up¬ 
per bound derived from the knowledge of the SR characterizing 
the AO correction of the atmospheric blur during the observa¬ 
tion. Moreover, non-negativity and normalization provide ad¬ 
ditional constraints. In conclusion, the nonconvex optimization 
problem we are considering can be formulated as follows 


min Jo(f,Ki,...,Kp;g,b) 
s.t. / > 0 , ^/(n) = c ; 

nG5 


(3) 


0<K,<5,-,yKy(n)=l; 7=l,...,j 


where sj is the upper bound on the PSP Kj derived from the 
knowledge of the SR characterizing the acquisition of gj. In 
conclusion the data of the problem are {g, b) and s - 

Another important consttaint can be provided by the require- 
m ent of band-lim i ting o f the PSPs, which is used for instance 


Desidera et al. ( 20061) . Indeed the band of the PSP, i.e the 


in 

set in Pourier space where the Pourier transform of the PSP is 
not zero, is known and consists, in general, of a disc in the case 
of a single aperture telescope and a union of three discs in the 
case of a Pizeau interferometer with LBT, the central one being 
the band of the mirrors of LBT and the side-ones the r eplica s 


due to interferometry (see, for instance, Bertero et al. . 201 ih . 


This constraint on the Pourier transform of the PSP, together 
with normalization and upper and lower bounds, defines a con¬ 
vex set. Unfortunately, since we use projection methods, the 
projection of an array on this set is not easily computable, even 
if methods for computing the projection on the intersection of 
convex sets are available. The difficulty is that these methods 
are not efficient and therefore can lead to an excessive compu¬ 
tational cost since the projection should be computed several 
times in the used iterative m ethods. However, f rom the numeri¬ 
cal experiments described in Prato et d] (l2013h we deduce that 
a suitable initialization of our algorithm, based on a PSP sat¬ 
isfying the band constraint, may lead to a set of reconstructed 
PSPs whose bands are very close to the desired ones. 

2.2. Blind deconvolution: alternating minimization 

Although the previous formulation of blind deconvolution re¬ 
quires the minimization of a nonconvex objective function, the 
constraints have a nice separable structure, since they involve 
separately the blocks of variables, defining a feasible convex 
set for each of them. In addition, the function is convex with 
respect to the different blocks of variables. In these settings, the 
solution of problem (O can be sought by means of an alternat¬ 
ing minimization (AM) strategy. 

The basic idea of AM is the cyclic minimization for the 
constrained problem with respect to one block of variables, 
updating its value for the next minimization step. This 
kind of approach is known in the literature also as nonlin¬ 
ear Gauss-Seidel or block coordinate descent method and its 
theoretic al propert i es have been deeply studied in the last 
decades ( Bertsek ai, 1 9991 ; Grippo and Sciandrone , 19991 2000l ; 


Luo and Tsenglll992h . 

In our case, each iteration of the AM method consists in solv¬ 
ing the following p \ constrained minimization problems of 
convex functions 


/ 


- arg min^,n-/o(/, Kf\ g, b) 

ik). 


(4) 


K\ = arg min^gjj 




]^{M) _ 

P 


= arg minji-eo ..., K, g, b) , 


neS 


where k is an index running on the AM iterations, Q is the set 
of the constraints on the object / and Q.j is the set of the con¬ 
straints on the PSP Kj. The limit points of the sequence gener¬ 
ated by this iteration scheme are also stationary points for the 
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constrained problem if each parti al problem is the minimiza¬ 
tion of a strictly convex function ( Bertsekas . 1999ll . This con¬ 
dition is not only sufficient but also necessary when more than 
two blocks of variables are involved. Indeed, [Powelll (1197311 


showed a counterexample where the strict convexity is not satis¬ 
fied, three blocks of variables are involved and the AM method 
fails to locate stationary points. 

Even when the hypothesis of strict convexity on each block 
of variables holds true, the convergence of the AM scheme can 
be proved only if each partial minimization problem is solved 
exactly, which is often impractical or too costly (in general it¬ 
erative methods are used). Quite surprisingly the convergence 
result is obtained without the assumption of strict convexity if 
SGP is used for solving inexactly the pa rtial minimizatio n prob¬ 
lems. This important result is proved in Bonettini ( 201 ih and is 
basic for the proposed approach to blind deconvolution. 

In conclusion, at each line of Eq. dUl the minimization is 
replaced by a given number of SGP iterations. These will be 
called inner iterations while the iterations of the AM scheme 
will be called outer iterations. Therefore the sequence gener¬ 
ated by the method depends both on the initialization of the 
first outer iteration (in the subsequent outer iterations the inner 
iterations are initialized with the results derived from the previ¬ 
ous one) and on the given numbers of inner iterations for each 
block of variables. In the present application it seems quite nat¬ 
ural to choose the same number of inner iterations for all the 
PSP’s blocks. 

2.3. SGP algorithm 

Since it is basic for the solution of the partial minimization 
problems it may be useful to briefly recal l the main points of 


the SGP algorithm (iBonettini et all l2009h even if its applica 


tion to astronomical imaging has already been described else 


where ( Prato et al. . 2012 : Bonettini and Prato . 2010l 20141) . To 
this purpose we remark that each minimization problem in the 
iterative scheme of Eq. (|4li has the following structure 


minyo(^) , 

heQ. 


(5) 


where, for simplicity, we omitted the dependence on the other 
variables and Q is the closed and convex set defined by the con- 
straints. The main difference with respect to Prato et aTI ( 2012 ) 
is that Q is a subset of the non-negative orthant defined by a 
suitable equality constraint. Therefore the projection on this set 
is more complex than that on the non-negative orthant. 

The main step of SGP is the computation of the k-th feasi¬ 
ble descent direction (where k is an index running on the inner 
iterations of a given AM iteration) 

- ar£>^V7o(/i®)) - /i® 

k 

by performing the following steps: 

a) The direction provided by the negative gradient 

is modified by a diagonal scaling matrix Dh with positive 
entries, which in all the subproblems of one AM iterations 
is given by 


Dk = diag(min(L2,max(Li,li®)) 


(Li, L 2 ) being given constants estimated from the extreme 
values of the image. 

b) A point on the scaled gradient direction is selected by 
choosing a multiplicative factor by means of an al- 
ternation of the generalized Barzilai-Borwein (BB ) rules 
( Barzilai and Borwein . 1988 : Bonettini et al. . 2009h 


JBBI) 


a, 


.(BB2) _ 




(7) 


where = /i® - and = V/q)^®) - 

and a suitable introduction of upper and lower 

bounds. 

c) The resulting point is brought back in the feasible set Q. 
by means of the projection Pq jj-i associated to the norm 
induced by D,', i.e. 


Cl.DT 


{h)^&xgmmy^n{h-yYD^\h-y). ( 8 ) 


Since the feasible sets of both the object and the PSPs in¬ 
volve a given number of inequalities plus an equality con¬ 
straint, in all cases we used a se cant-based routine devel¬ 
oped by Dai and Pletcher ( 2006h . which is able to compute 
the projection with a computational cost gro wing linearly 


in tim e with respect to the image size (see also lPrato et al 
2013h. 


2.4. Discussion 

In the previous approach, the blind deconvolution problem 
is formulated as the constrained minimization of a nonconvex 
function which depends on an extremely large number of vari¬ 
ables, about 10^ in the numerical experiments described in this 
paper. Since the constraints used in our approach imply that the 
sequences of objects and PSPs generated by the inexact AM 
method are bounded, it follows that these sequences have limit 
points. We can add that, even if it is difficult to provide a the¬ 
oretical evidence of the existence of a unique limit point, in all 
our numerical experiments the sequences produced by the in¬ 
exact AM method have a convergent behaviour. Ho wever, ac¬ 
cordin g to the general convergence result proved in IBonettini 
JMH) , we can only state that the limit points are stationary 
points of the function, hence not necessarily minimizers. 

As far as we know, there is no practical way for establishing 
if these points are minimizers or not. In fact, it should be nec¬ 
essary to manage the Hessian of the function in these points but 
this is an absolutely intracta ble matrix even if one can write it 


explicitly ( Prato et al. . 20131) . Since one can use different ini- 


(6) 


tializations of the iterative procedure and different numbers of 
inner iterations and these different choices can produce differ¬ 
ent results, in a practical application we do not see an approach 
better than that of doing different attempts and look for that 
providing the most sensible solution. 
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An additional difficulty is that it may happen, as we show 
by some numerical simulations, that a sensible solution cor¬ 
responds to a value of the objective function which is greater 
than the value of the same function corresponding to a solution 
which is clearly unphysical. It is obvious that these situations 
should not be surprising because the problem of blind decon¬ 
volution is nonconvex and therefore the objective function can 
have several local minimizers as well as stationary points. Since 
the objective function has a simple structure it should be impor¬ 
tant to characterize the sets of these points, but an approach to 
this problem presently is not available, as far as we know. 

The advantage of our method is that it is mathematically 
sound, it provides sequences with limit points, very frequently 
with a unique limit point and therefore, if the user is conscious 
of the difficulties of the problem, he can attempt to use this 
method for obtaining different solutions in practical applica¬ 
tions and select that looking as the most appealing. In the next 
section we attempt to provide a few hints for helping the user 
in the choice of the parameters of the method and, in particular, 
of its initialization. 


3. Image simulation 

We model the im ages according to the model proposed in 
Snvder et al.l (119941) for images acquired with a CCD camera, 
i.e. each pixel is affected by background (due to sky emission, 
dark current, etc.), photon counting noise (described by a Pois¬ 
son distribution) and additive read-out noise (RON) described 
by a Gaussian distribution. 

If the RON variance is cr^, in the deconvolution process it 
can be approximated by a Poisson distribution with parameter 
cr^ if cr^ is added both to the detected ima ges and the corre¬ 


sponding backgrounds (ISnvder et al.LIl995h . Therefore all the 


pixel values of the detected images can be viewed as realiza¬ 
tions of suitable Poisson random variables if in Eq. (ID we in¬ 
tend that gj, bj have been modified according to this approach. 
Therefore, in our numerical simulations we perturb the images 
with Poisson and additive Gaussian noise but in the deconvo¬ 
lution algorithms we use the images and backgrounds modified 
as above. 

All the images and the PSFs considered in our numerical ex¬ 
periments are sized 256 x 256 pixels in the single image case, 
with a pixel size of 15 mas, and 512 x 512 pixels in the multiple 
image case, with a pixel size of 5 mas. Moreover all images, 
except one indicated in Sect. 4.1.1, are obtained by adding 10 
frames in order to avoid saturation of the detector, as we discuss 
in the following, so that the variance of the RON will be 10 cr^. 


3.1. Single image simulation 

In this case we use two PSFs in K-band with SR = 0.81 and 
0.62 respectively, modeling the optics of a single mirror of LET, 
with diameter 8.4 m, and the effect of the adaptive optics system 
FLAG using the power spectrum of the wavef ront residual of 


the A O correction as measured at the telescope (lEsposito et al. 


20121) . To the noise-free image, obtained by convolving the ob¬ 
ject with one of these PSFs, a background in K-band is added 


and the result is corrupted with Poisson and additive Gaussian 
noise. In order to avoid saturation of the detector (a maximum 
number of 5 x 10^ photons per pixel is assumed in a single 
frame) the image is obtained by co-adding n frames. More pre¬ 
cisely, in the case of a stellar system the procedure for image 
generation is the following. 

• We establish the coordinates of the stars and we fix their 
magnitudes in K-band. 

• We compute the integration time which does not produce 
saturation of the detector by taking into account the col¬ 
lection area of the telescope, the overall efficiency of the 
acquisition system (assumed equal to 30%), and the flux of 
the brightest star multiplied by the peak value of the PSF. 
This is the integration time of a single frame and is used 
for computing the number of frames n required for obtain¬ 
ing an acceptable SNR for all the stars of the system. 

• We generate noise-free images by shifting, with sub-pixel 
precision, the PSF to the positions of the stars and adding 
these shifted PSFs, each one weighted with a weight corre¬ 
sponding to the magnitude and the total observation time. 

• These images are perturbed by adding a background in K- 
band, corresponding to about 13.5 mag arcsec"^, and by 
corrupting the results with Poisson and additive Gaussian 
noise (RON); the variance of the RON is ncr^, thus corre¬ 
sponding to the RON of n frames; we take cr = 10 e / px. 

3.2. Multiple image simulation 

As concerns the simulation of LN images, we recall that the 
instrument combines in a Fizeau mode the beams coming from 
the two mirrors of LET whose center-to-center distance is about 
14.4 m. Therefore the maximum baseline available is 22.8 m 
and the resolution achievable by a single LN image is that of a 
22.8 m telescope in the direction of the baseline and that of a 8.4 
m telescope in the orthogonal direction. For a given orientation 
the PSF of LN looks as that of a 8.4 m telescope, modulated by 
the interference fringes, orthogonal to the direction of the base¬ 
line. In order to get a more uniform resolution one must acquire 
and combine different images with different orientations of the 
baseline. 

It is important to remark that the orientation of the fringes 
does not depend on the orientation of the baseline because the 
camera is rotating with the baseline and therefore the fringes 
have always the same direction (for instance the vertical one) 
in the image array. In other words two images of the same 
scientific object with two different orientations of the baseline 
correspond to two rotated versions of that object. This specific 
feature implies that one should introduce rotation matrices in 
the formulation of the problem. However we verified that the 
computation of hundreds or thousands of rotations in hundreds 
or thousands of inner iterations introduces large computational 
errors. Therefore we considered the approach which consists 
in derotating the images in such a way that they correspond 
to aligned versions of the object /. The price to be payed is 
that the derotation of discrete images modifies their statistical 
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properties. In order to estimate this effect we considered the 
rotation of a constant array perturbed by Poisson noise. We 
found the following results: 

• before rotation the histogram of the array is a Gaussian 
with the same mean and variance; after a rotation based on 
spline interpolation the histogram is still a Gaussian with 
the correct mean but a smaller variance; 

• the support of the autocorrelation of the rotated image is a 
3x3 square; 

• if we use a different rotation approach which consists in 
attributing the value of a pixel before rotation to the pixel 
with maximum overlapping after rotation (nearest neigh¬ 
bor approximation), the statistics is preserved but the qual¬ 
ity of the image is degraded. 


As a consequence of this analysis we decided to use in the ap¬ 
proach derotated images. 

The procedure adopted in our numerical experiments is sim¬ 
ilar to that used in the case of a single image. We consider 
two sets of PSFs in K-band with SR respectively 0.77 and 0.46, 
corresponding to orientation angles of the baseline indicated as 
0°, 60° and 120°, all with vertical fringes (for simplicity we 
take the same SR for the three orientations). The hrst PSF of 
each se t has been generated by m eans of the software package 


LOST (lArcidiacono et al.L 120041) . the second by reflecting the 


first one with respect to the central line and the third by taking 
the arithmetic mean of the first two. In this way the three PSF of 
each set have exactly the same SR. Then the generation of the 
corresponding LN images is similar to that of the single image 
case by modifying the first item as follows. 


• We establish the coordinates of the stars corresponding to 
the observation at 0° and we compute, with sub-pixel pre¬ 
cision, their coordinates if the system is rotated by 60° and 
120° respectively. 


The rest of the procedure is unchanged and applied to the three 
images but at the end we must add the following item. 

• The images corresponding to 60° and 120° are derotated 
in order to align the object in the three images and three 
arrays containing the object are extracted from the full im¬ 
ages. 

The derotated images are used in the definition of the objec¬ 
tive function and in the blind algorithm, which therefore will 
produce derotated PSFs. 


the case of a stellar system we consider a magnitude average 
relative error (MARE) defined by 

1 |m; - m,j 

MARE = - > !-4 ^ (9) 

q Mi 

where q is the number of stars and m,, to, are respectively the 
reconstructed and the true magnitudes. 

As concerns PSF reconstruction, in the case of single im¬ 
age we consider the root-mean-square error with respect to the 
true one, defined as usual in terms of the {2 norm of their dif¬ 
ference. In the case of LN images generated according to the 
previous procedure, since the blind algorithm produces a set 
of three PSFs, two of them being derotated with respect to the 
ones used for generating the images, for comparison we must 
derotate the original ones. If we denote as Kj the derotated orig¬ 
inal PSF, then we measure the quality of the reconstruction by 
means of the root-mean-square error (RMSE) 


\\Kj-Kj\\ 

WKjW 


( 10 ) 


where Kj is the reconstructed PSE and || ■ || denotes the usual 
f’ 2 -norm. 


4.1. Binary systems 

We first consider the simple case of binary systems. More 
precisely we consider nine cases by varying both separation 
and magnitude of the stars. By keeping fixed the magnitude 
of the primary, i.e. nii = 15, we take for the magnitude of the 
secondary m 2 - 15, 16 and 17. Moreover for each choice we 
consider three possible angular separations: d - 60, 120 and 
240 mas in the single image case and d = 20, 40 and 80 mas 
in the LN case. In both cases the first separation corresponds 
to the resolution limit of the instrument while the last is four 
times larger. In all cases, as described in the previous section, 
we compute the integration time of a frame in such a way that 
the number of counts in the image pixel corresponding to the 
position of the primary does not exceed 5x10^. As stated in 
the previous section, we consider 10 frames per image, both in 
the single and in the multiple image case, so that the peak value 
of the photons is about 5 x 10^ for all images. Since in the case 
of LN we have three images, in this case the SNR is higher than 
in the single image case. 

In Eig. 1 we show the images of the binaries with mi - m 2 - 
15 and different angular separations; in the hrst row those of the 
single image case and in the second row those of the multiple 
image case corresponding to the 0° baseline, all obtained with 
the PSE with the highest SR. The difficulty in reconstructing the 
binary with separation d -2Q mas is obvious. 


4. Numerical results 

In order to evaluate the quality of the reconstructions ob¬ 
tained with our blind method we need some hgures of merit. 

As concerns the reconstruction of a binary we consider the 
relative absolute error on the magnitudes of both stars while in 


4.1.1. Single image 

Eor the convenience of the reader we give the computed inte¬ 
gration time avoiding saturation in a single frame: 40 sec for SR 
= 0.81 and 52 sec for SR = 0.62. As already stated the images 
are obtained by adding 10 frames. These are the input images of 
the blind algorithm together with the value of the background. 
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Figure 1: Examples of input images of binaries with magnitudes m[ = m 2 = 15. In the first row those of the single image case, con'esponding to PSF with SR = 
0.81: from left to right, angular separation of 60, 120 and 240 mas. In the second row those of the multiple image case, con'esponding to the PSF with SR = 0.77: 
from left to right, angular separation of 20, 40 and 80 mas. These images correspond to the first orientation of the baseline and only the central part of the images 
256 X 256 is displayed. In the two other orientations the binaries appear rotated by 60 and 120 degrees respectively. Images are displayed in log scale. The length 
coiTesponding to 0.5 ai'csec is also indicated. 


d = 60 mas ; m2=15 



d=60 mas ; m2=15 



Figure 2: Behaviour, as a function of the number of iterations, of the normalized objective function (left panel) and of the RMSE on the PSF (right panel). The 
pai'ameters of the binary are indicated in the figure. The plots refer to the PSF with SR = 0.81. 
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In a firs t attem pt we use the initialization already used in 
Prato et aT] ( 2013 ) and in other papers, namely a constant array 


for the object and the autocorrelation of the diffraction-limited 
PSF for the PSF. Indeed, this initialization has produced very 
promising results in our previous paper, where a much higher 
SNR was assumed. We use 1000 outer iterations in the case SR 
= 0.81 and 2000 outer iterations in the case SR = 0.62. Indeed 
in the case of a lower SR we have a lower quality of the images 
and, presumably, a larger number of iterations is re quired. As 
concerns the inner iterations, as in Prato et al.l ( 2013 ) we use 50 
SGP iterations for the object and one SGP iteration for the PSF. 


In Table 1 we give the results obtained with the previous 
choice. As a first remark, the binaries and the PSFs are re¬ 
constructed satisfactorily in all cases except the closest binaries 
id - 60 mas) with different magnitudes. Indeed, the indica¬ 
tion 100% in the column for Am 2 lm 2 means that the method 
reconstructs only one star, which sometimes is not exactly in 
the position of the primary but slightly shifted in the direction 
of the secondary. Since its magnitude is computed using a 3x3 
square centered on the true position of the primary, the error on 
its magnitude is, in general, not too large. On the other hand the 
error on the PSFs is very large, as one should expect since the 
secondary is missed. This point deserves further investigation. 


In Fig. 2 we show, in a particular case, the behaviour of the 
normalized objective function, defined by UqIN^ with Jq given 
in Eq. ([TJ (with p - 1), and of the RMSE on the PSE as func¬ 
tions of the number of iterations. Similar behaviors are obtained 
in all cases where a sensible result is obtained. This result sug¬ 
gests that presumably convergence is reached after 1000 itera¬ 
tions even if, as previously discussed, it is difficult to establish 
numerically the convergence of a sequence. 


A second remark is that, according to statistical properties of 
Poisson random variables, if we compute the value of the nor¬ 
malized objective function by inserting in Eq. ([T]i the noisy and 
th e noise-free images we should obtain a value very close to 
1 ( Bertero et al.L 2010; Zanella et al.L 2009ll . This is just what 
we obtain using our simulated images (this result also demon¬ 
strates the accuracy of the approximation of the RON with a 
Poisson random variable). However the limiting values of the 
normalized objective function obtained in our experiments are 
definitely smaller than 1, an effect already remarked in our pre¬ 
vious paper. 


Coming back to the problem of the unresolved binaries, we 
point out that, if we deconvolve the images using the PSP used 
for their generation {inverse crime) all the binaries are correctly 
reconstructed with small errors on their magnitudes. Therefore 
the failure of our experiment may be due to a failure of the 
method or to an inappropriate initialization or to inappropriate 
choices of the internal iterations. 


Several attempts with different numbers of internal iterations 
did not improve the results. Therefore we searched for an initial 
PSF with a SR value closer to the correct one and with the prop¬ 
erty of being band-limited with the band of the LET mirror. A 
possible choice is obtained by means of the diffraction-limited 
PSF of LET, let us say K, by looking for an initial guess /T® of 


the following form 


1 + cj 

which is band-limited and satisfies the normalization condition. 
The constant a> should be selected in such a way that has 
the correct SR value, i.e. max - SR max (K). We obtain 

(SR max{K) - l)w = (SR - 1) max(.^) (12) 

and, by neglecting 1 with respect to the first term in the l.h.s. of 
this equation, we obtain cu = (1 - SR)/(SR N^). 

The results obtained with this initialization, using again 50 
SGP iterations for the object and one for the PSF, are reported 
in Table 2. Since the convergence is slower than in the previ¬ 
ous case we use 2000 outer iterations for SR = 0.81 and 3000 
iterations for SR = 0.62. 

Ey comparing the results reported in the two tables we re¬ 
mark that the two different initializations provide very similar 
results in all cases where they succeed or they fail; in other 
words they provide sequences of iterations which presumably 
converge, even if with a different rate, to the same point, which 
is a stationary point of the objective function. Obviously we 
believe that it is also a minimizer. In the case of separation 60 
mas and m 2 - 16 the algorithm, equipped with the new ini¬ 
tialization, is able to reconstruct the binary and the PSF with a 
satisfactory accuracy for both values of SR. We remark that the 
value of the objective function is higher than that correspond¬ 
ing to the result provided by the first initialization, which is not 
correct. This fact clearly indicates the existence of several sta¬ 
tionary points or minimizers or both. Of course it should be 
nice to establish that the result of the first initialization is a sta¬ 
tionary point and that of the second a minimizer; but, as already 
remarked such a verification is practically impossible. Finally, 
in the case 1112 = 17 also the new initialization is unable to pro¬ 
vide the correct results. 

The results obtained in the multiple image case and described 
in the next subsection suggest that this negative result may be 
due to an insufficient value of the SNR. Therefore, in the case 
m 2 = 17 we generated an image which is the sum of 30 frames 
(we point out that, as already remarked, in the considered mul¬ 
tiple image case we have three times the photons of the single 
image case). Using again 2000 iterations, we find that the al¬ 
gorithm, with the second initialization, can resolve the binary 
in the case SR = 0.81 (even if with a large reconstruction error, 
about 9 %, on the PSF) but not in the case SR = 0.62. 

However in these difficult cases we observe a new phe¬ 
nomenon; even if in the limit the results are not satisfactory, 
the PSF reconstruction error exhibits a minimum before con¬ 
vergence. If we consider the reconstructions corresponding to 
these minima, then, in the case of the first initialization, the 
minima do not correspond to a situation where the binary is 
resolved. On the other hand, in the case of the second initial¬ 
ization, the binary is resolved for both SR values, with a 2.03 
% PSF error in the case SR = 0.81 (574 iterations) and a 7.13 
% error in the case SR = 0.62 (1739 iterations). Such a result 
presumably indicates the need of introducing a regularization 
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Table 1: Single image case - Binary reconstructions provided by the algorithm initialized with the autocorrelation of the dilfraction-limited PSF. In the first column 
the value of the SR, in the second the angular separation, in the third the magnitude of the secondaiy, in the fourth and fifth the errors on the magnitudes of the two 
stars. In the subsequent column we give the RMSE for the reconstructed PSF. Finally in the last two columns we give the value of the normalized objective function, 
defined by UqIN'^, as computed at the end of the iterations, and the number of outer iterations. 


SR 

d (mas) 

1712 

Ami 1 mi 

Am2/m2 

RMSE 

Jq norm. 

IT 



15 

0.05% 

0.03% 

0.94% 

0.6071 

1000 


60 

16 

2.37% 

100% 

40.41% 

0.5549 

1000 



17 

0.99% 

100% 

16.77% 

0.5964 

1000 



15 

<0.01% 

<0.01% 

0.76% 

0.6047 

1000 

0.81 

120 

16 

0.03% 

<0.01% 

1.11% 

0.6296 

1000 



17 

0.03% 

0.17% 

1.40% 

0.6254 

1000 



15 

<0.01% 

<0.01% 

0.79% 

0.5999 

1000 


240 

16 

0.02% 

0.03% 

0.83% 

0.6273 

1000 



17 

<0.01% 

0.04% 

1.17% 

0.6229 

1000 



15 

0.18% 

<0.01% 

1.08% 

0.5338 

2000 


60 

16 

2.31% 

100% 

34.37% 

0.4635 

2000 



17 

1.05% 

100% 

16.87% 

0.4983 

2000 



15 

0.15% 

0.14% 

1.04% 

0.5261 

2000 

0.62 

120 

16 

0.02% 

0.01% 

1.28% 

0.5419 

2000 



17 

0.04% 

0.25% 

1.59% 

0.5329 

2000 



15 

0.04% 

0.04% 

1.00% 

0.5309 

2000 


240 

16 

<0.01% 

0.06% 

1.13% 

0.5537 

2000 



17 

0.05% 

0.36% 

1.80% 

0.5361 

2000 


Table 2: Single image case - Binary reconstructions provided by the algorithm initialized with the diffraction-limited PSF plus a constant selected for satisfying the 
SR constraint (see the text). The structure of the Table is the same of Tabled 


SR 

d (mas) 

1712 

Ami 1 mi 

Am2/m2 

RMSE 

Jq norm. 

IT 



15 

0.02% 

0.04% 

0.82% 

0.6067 

2000 


60 

16 

0.09% 

0.16% 

2.05% 

0.6237 

2000 



17 

1.01% 

100% 

17.08% 

0.5956 

2000 



15 

<0.01% 

<0.01% 

0.77% 

0.6046 

2000 

0.81 

120 

16 

0.02% 

<0.01% 

1.09% 

0.6294 

2000 



17 

0.02% 

0.15% 

1.35% 

0.6253 

2000 



15 

<0.01% 

<0.01% 

0.80% 

0.5989 

2000 


240 

16 

0.02% 

0.02% 

0.82% 

0.6271 

2000 



17 

<0.01% 

0.02% 

1.12% 

0.6227 

2000 



15 

0.02% 

<0.01% 

1.11% 

0.5333 

3000 


60 

16 

0.12% 

0.25% 

2.64% 

0.5354 

3000 



17 

1.05% 

100% 

16.87% 

0.4983 

3000 



15 

0.01% 

0.01% 

1.06% 

0.5258 

3000 

0.62 

120 

16 

0.02% 

<0.01% 

1.26% 

0.5419 

3000 



17 

0.04% 

0.25% 

1.58% 

0.5329 

3000 



15 

0.03% 

0.03% 

0.99% 

0.5304 

3000 


240 

16 

<0.01% 

0.06% 

1.12% 

0.5537 

3000 



17 

0.05% 

0.36% 

1.80% 

0.5361 

3000 
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Figure 3: Single image case - PSF reconstruction in the case of the binary with d = 60 mas and m 2 = 17. The input image is the sum of 30 frames (see text). 
These PSFs correspond to the minima of the reconstruction error. First column: the true PSF with SR = 0.81 (top) and SR = 0.62 (bottom). Second column: 
PSF reconstruction provided by the algorithm initialized with the autocorrelation of the diffraction-limited PSF. Last column: PSF reconstruction provided by the 
algorithm initialized with the diffraction-limited PSF plus a constant. In each panel we also show a zoom of the core of the PSF which makes evident artifacts due 
to the secondary. All images are displayed in log scale. 


of the PSF in the objective function, at least for treating the 
most difficult cases. In Fig. 3 we show the reconstructions of 
the PSF corresponding to the minimum reconstruction errors. 
Artifacts due to the missed secondary are visible in the case of 
the first initialization and also in the case SR = 0.62, since the 
reconstructed secondary is fainter than the true one. 

4.1.2. Multiple images 

In this case the integration time of a nonsaturated frame is 95 
sec for SR = 0.77 and 167 sec for SR = 0.46. For each binary 
and orientation angle we consider again 10 frames, so that we 
have approximately the same number of photons in all images. 

We preliminarily remark that, if we compute the value of the 
normalized objective function (which is now given by 2JqI3N^) 
by inserting in Eq. O the noisy and the noise-free images be¬ 
fore derotation, we expect to obtain a value very close to 1 and 
this is just what we obtain. But this is not true if we compute 
the same quantity using the derotated images. Indeed, for the 
nine binaries as well as for the other objects, we always obtain 
a smaller value, namely 0.63. Since this value is independent of 
the object and PSFs, this effect is clearly due to the modification 
of the statistical properties of the data introduced by the derota¬ 
tion, as briefly discussed in Sect. 3.2. In any case the limiting 
values of the normalized objective function obtained in our ex¬ 
periments are definitely smaller than the values corresponding 
to the input objects and images, an effect already remarked in 
the previous case. 

As in the single image case we first use as initialization a con¬ 
stant array for the object and the autocorrelations of the ideal 
PSFs for the three PSFs. The results of the reconstructions ob¬ 
tained with this initialization are reported in Table 3. We obtain 
that only when both stars have the same magnitude the method 
is able to reconstruct both the binary and the PSFs with suf¬ 


ficient accuracy. When we have different magnitudes for the 
two stars the method is in general failing to reproduce the sec¬ 
ondary, except in the case of separation d - SO mas; in this 
case a binary with difference of magnitude Am = 1 is also re¬ 
constructed. As in the single image case, the indication 100% 
in the column for Am 2 /m 2 means that the method reconstructs 
an object which contains only one bright star (in one case the 
centroid is shifted one pixel in the direction of the secondary. 
These results show that, even if we have a higher SNR as al¬ 
ready discussed, the multiple image case is more difficult than 
the single one. 


If we deconvolve the derotated images using the derotated 
PSFs (this is not exactly an inverse crime because the images 
were generated with non derotated PSFs) all the binaries are 
correctly reconstructed with small errors on the magnitudes. 
Therefore the failure of our experiment may be due again to 
an inappropriate initialization (the autocorrelations of the ideal 
PSFs have a SR value of about 0.35, much smaller than the 
SR of the PSFs used in image generation) or to inappropri- 
ate choices of the internal iterations. Also in this case, as in 


Prato et al.l (120131) and in the single image case, we use 50 SGP 


iterations for the object and one SGP iteration for each PSF. 
However several attempts with different numbers of internal it¬ 
erations did not improve the results. Therefore, as in the single 
image case, we use as a new initialization of the PSFs the ideal 
PSFs of LN with the addition of a small constant selected in 
such a way to satisfy normalization and SR value. The results 
obtained with this initialization, using again 50 SGP iterations 
for the object and one for the PSFs, are reported in Table 4. 
Since the convergence is slower than in the previous case we 
use 2000 outer iterations. 


With the new initialization the blind method succeeds in re- 
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Table 3: Multiple image case - Binary reconstructions provided by the algorithm initialized with the autocorrelations of the ideal PSFs. In the first column the value 
of the SR, in the second the angulai* separation, in the third the magnitude of the secondaiy, in the fourth and fifth the eiTors on the magnitudes of the primary and 
the secondary star. In the subsequent three columns we give the RMSE for the three PSFs. Finally in the last two columns we give the value of the normalized 
objective function, defined by 27o/3A^^, as computed at the end of the iterations, and the number of outer iterations. 


SR 

d (mas) 

m2 

Am\lmi 

AOT 2 /OT 2 

RMSEo» 

RMSE6o° 

RMSE 120 " 

Jo norm. 

IT 



15 

0.28% 

0.23% 

1.52% 

2.50% 

1.87% 

0.2241 

1000 


20 

16 

2.12% 

100% 

30.00% 

31.14% 

21.79% 

0.1706 

1000 



17 

0.87% 

100% 

14.93% 

15.38% 

11.00% 

0.1890 

1000 



15 

0.30% 

0.28% 

1.66% 

1.84% 

2.56% 

0.2647 

1000 

0.77 

40 

16 

2.20% 

100% 

41.65% 

33.98% 

33.80% 

0.2049 

1000 



17 

0.58% 

100% 

14.78% 

14.48% 

15.90% 

0.1954 

1000 



15 

0.20% 

0.21% 

1.09% 

0.83% 

0.83% 

0.2180 

1000 


80 

16 

0.18% 

0.23% 

1.27% 

0.96% 

0.99% 

0.2164 

1000 



17 

0.87% 

100% 

19.28% 

19.06% 

19.04% 

0.1893 

1000 



15 

1.27% 

1.19% 

9.67% 

9.69% 

11.39% 

0.1335 

1000 


20 

16 

0.29% 

100% 

32.61% 

33.10% 

29.37% 

0.0836 

1000 



17 

0.18% 

100% 

13.91% 

14.12% 

12.54% 

0.0795 

1000 



15 

0.89% 

0.88% 

5.15% 

5.62% 

5.64% 

0.1516 

1000 

0.46 

40 

16 

0.27% 

100% 

47.22% 

41.54% 

36.75% 

0.1360 

1000 



17 

0.33% 

100% 

13.91% 

14.67% 

14.31% 

0.0944 

1000 



15 

0.68% 

0.68% 

3.05% 

2.60% 

2.58% 

0.1042 

1000 


80 

16 

0.52% 

0.60% 

1.87% 

1.43% 

1.44% 

0.0850 

1000 



17 

0.55% 

100% 

15.99% 

15.97% 

15.98% 

0.0883 

1000 


Table 4: Multiple image case - Binary reconstructions provided by the algorithm initialized with the ideal PSFs plus a constant selected for satisfying the SR 
constraint (see the text). The structure of the Table is the same of Tabled 


SR 

d (mas) 

m 2 

Ami/mi 

AOT 2 /OT 2 

RMSEo» 

RMSE6o° 

RMSE 120 " 

Jo norm. 

IT 



15 

0.44% 

0.34% 

2.86% 

4.23% 

3.42% 

0.2277 

2000 


20 

16 

0.27% 

0.21% 

1.56% 

1.82% 

1.74% 

0.2209 

2000 



17 

0.07% 

1.10% 

2.53% 

2.70% 

1.78% 

0.2095 

2000 



15 

0.45% 

1.03% 

5.47% 

4.61% 

6.73% 

0.2670 

2000 

0.77 

40 

16 

0.25% 

0.39% 

1.63% 

2.85% 

2.76% 

0.2220 

2000 



17 

0.11% 

0.73% 

2.05% 

2.78% 

2.86% 

0.2102 

2000 



15 

0.35% 

0.35% 

2.28% 

1.51% 

2.23% 

0.2204 

2000 


80 

16 

0.25% 

0.26% 

1.32% 

1.06% 

1.14% 

0.2179 

2000 



17 

0.19% 

0.40% 

1.32% 

1.02% 

0.99% 

0.2125 

2000 



15 

0.80% 

0.57% 

4.02% 

4.51% 

6.83% 

0.1037 

2000 


20 

16 

0.38% 

0.95% 

3.76% 

3.95% 

2.52% 

0.0811 

2000 



17 

0.07% 

6.32% 

9.05% 

10.33% 

6.29% 

0.0697 

2000 



15 

0.64% 

2.18% 

11.23% 

6.72% 

8.73% 

0.1409 

2000 

0.46 

40 

16 

0.49% 

0.81% 

2.21% 

3.20% 

2.99% 

0.0837 

2000 



17 

0.02% 

5.89% 

8.70% 

7.68% 

8.84% 

0.0716 

2000 



15 

0.56% 

0.55% 

8.54% 

4.66% 

4.64% 

0.1100 

2000 


80 

16 

0.57% 

0.53% 

1.99% 

1.48% 

1.49% 

0.0846 

2000 



17 

0.48% 

0.92% 

2.36% 

1.95% 

1.96% 

0.0785 

2000 
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Figure 4: Multiple image case - PSF reconstruction in the case of the binary with d = SO mas and m 2 = 17. First column: the true PSF with SR = 0.77 (top) 
and SR = 0.46 (bottom). Second column: PSF reconstruction provided by the algorithm initialized with the autocorrelations of the ideal PSFs. Last column: PSF 
reconstmction if the algorithm is initialized with the ideal PSFs plus a constant. All images (only the central part 256x256 is shown) are displayed in log scale and 
coiTespond to the first orientation of the baseline. 


constructing all the binaries with sufficient accuracy as well as 
the PSFs. We can add that in most cases both the normalized 
objective function and the RMSE on the PSFs have a convergent 
behaviour while, in a few cases, the errors are still decreasing 
after 2000 iterations, thus indicating that a larger number of it¬ 
erations could still improve the solution. A comparison of the 
values of the objective function reported in the two tables shows 
that, in some of the cases where the first initialization is failing, 
the values in Table 3 are smaller than the corresponding val¬ 
ues in Table 4. This phenomenon was already observed in the 
single image case and means that different stationary points or 
minimizers are present. 

A few more comments on the two tables. If one looks care¬ 
fully at the reported results one can remark that, even if the re¬ 
sults obtained with the second initialization are globally better 
than those obtained with the first one, this may not be true for 
particular cases (compare, for instance, the results for d - 40 
mas and Am = 0). Moreover, the errors obtained with the sec¬ 
ond initialization do not vary in a regular way with the variation 
of angular distance and difference of magnitude. These behav¬ 
iors can be due to the fact that 2000 iterations may not be suffi¬ 
cient for assuring convergence of the method in the case of the 
second initialization. We did not push further the iterations be¬ 
cause in the case of three 512x512 images the computation time 
is considerable. By assuming possible fluctuations due to insuf¬ 
ficient number of iterations, a reasonable conclusion seems to 
be that, as in the single image case, the two initializations lead 
to the same limit point when the first one is successful. 

In Fig. 4 we show an example of reconstructions of the PSF 
at 0°, for both SR values, when the unknown object is a binary 
with = 80 mas and m 2 - 17. From the reconstructions dis¬ 
played in the second column and obtained by initializing with 
the autocorrelations of the ideal PSFs, it is evident that they 


contain a contribution coming from the secondary, while this 
contribution is practically absent in the reconstructions obtained 
with the other initialization and displayed in the third column. 


4.2. Star clusters 


In a second experiment we conside r two models of st ar clus¬ 
ter. The first is already considered in IPrato et alJ (120 13h and is 
based on an image of the brightest stars of the Pleiades open 
cluster; for this reason, we call it “open star cluster”. It consists 
of nine stars that we take, in this paper, with magnitudes rang¬ 
ing from 14.4 to 17.1. In the single image case, the minimum 
distance between two stars is 120 mas, while the maximum dis¬ 
tance is 1434 mas, with a mean distance of about 690 mas. In 
the multiple image case, considering the different pixel scale, 
we reduce of one third all the distances. 

As a second example we consider a model that we call “glob¬ 
ular star cluster”. For simplicity, only 150 stars are considered 
within the field of view, representing a very low crowding con¬ 
dition. The positions of the stars are randomly computed fol¬ 
lowing a Gaussian distribution around the center of the image 
(with a standard deviation of about 450 mas in the single image 
case and of about 150 mas in the multiple image case); similarly 
the magnitudes of the stars are randomly distributed around m 
= 16 with a standard deviation of about 0.4. It turns out that the 
brightest star of the cluster has m = 14.8. 

Again, we limit the maximum number of counts in each 
frame to 5 x 10"^, keeping fixed to 10 the number of frames. 
In Fig. 5 we show the images of the two star clusters provided 
by the PSFs with the highest SR. 


4.2.1. Single image 

In the case of the “open star cluster”, the integration time of 
a single frame is 22 sec for SR = 0.81 and 29 sec for SR = 0.62 
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Figure 5: Top panels: the input images of the “open star cluster” (left) and of the “globular star cluster” (right) in the single image case with SR = 0.81. Bottom 
panels: the input images of the two clusters in the case of SR = 0.77 and with 0° of the baseline (only the central part 256x256 is shown). All images are displayed 
in log scale. The length con'esponding to 0.5 arcsec is also indicated. 


Table 5: Single image case - The reconstruction errors in the case of the two star cluster models. In the first column the “open cluster” is labelled by OC, while 
the “globular cluster” is GC. In the second column, we give the value of SR, while in the third column we give the initialization of the algorithm, denoting by A 
the autocorrelation of the diflFraction-limited PSF and by C the diflFraction-limited PSF plus a constant selected for satisfying the SR constraint (see the text). In the 
subsequent columns, we give the value of the magnitude average reconstruction eiTor (MARE) defined in Eq. lE), and the RMSE for the reconstructed PSF. Finally 
in the last two columns we give the value of the normalized objective function, defined by 2JqIN^ , as computed at the end of the iterations, and the number of outer 
iterations. _ 


Star Cluster 

SR 

Init. 

MARE 

RMSE 

Jq norm. 

IT 


0.81 

A 

0.06% 

0.84% 

0.5993 

2000 

OC 

C 

0.06% 

0.85% 

0.5991 

5000 

0.62 

A 

0.09% 

1.22% 

0.5165 

4000 


C 

0.09% 

1.22% 

0.5165 

10000 


0.81 

A 

0.06% 

0.87% 

0.5123 

3000 

GC 

C 

0.06% 

0.82% 

0.4989 

5000 

0.62 

A 

0.07% 

1.07% 

0.4622 

6000 


C 

0.07% 

15.82% 

0.5698 

10000 
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while in the case of the “globular star cluster” these times are 
respectively 32 and 42 sec. 

We applied to the four images our blind algorithm using both 
initializations introduced in the case of the binaries. The results 
are reported in Table 5. In the case of the “open star cluster” 
and both values of SR the two initializations seem to provide 
sequences of iterations converging to the same point. If we 
look at the image shown in the upper left panel of Fig. 5 we 
can observe that it contains sufficiently well-separated star im¬ 
ages which can allow a good estimation of the PSF by the blind 
algorithm. 

The situation is a bit different in the case of the “globular star 
cluster” and we can understand this fact if we look at the upper 
right panel of Fig. 5. In the case of the higher SR value both 
initializations lead essentially to the same result. The small dif¬ 
ferences may be due to different convergence rates and could 
be removed by a more accurate tuning of the number of it¬ 
erations. On the other hand in the case of the lower SR ra¬ 
tio the first initialization, based on the autocorrelation of the 
diffraction-limited PSF, provides the best PSF reconstruction 
(also corresponding to a lower value of the objective function). 
It seems that the two initializations lead to two different station¬ 
ary points. In conclusion, for this particular object one can state 
that the first initialization may provide a better result than the 
second one. 

4.2.2. Multiple images 

In the case of the “open cluster” model, the integration time 
is 53 sec for SR = 0.77 and 93 sec for SR = 0.46. On the other 
hand the integration time for the “globular cluster” images is 78 
sec for SR = 0.77 and 136.5 sec for SR = 0.46. 

In both cases we apply our blind algorithm using the two ini¬ 
tializations already used in the previous sections, with 50 inner 
SGP iterations for the object and one iteration for each PSF. 
The results obtained for the “open cluster” with the two initial¬ 
izations are given in the first two rows of Table 6 in the case 
SR = 0.77 and in the following two rows those obtained in the 
case SR = 0.46. Similarly the results obtained for the “globular 
cluster” are given in the second half of the same table. 

In the multiple image case the situation is more complex than 
in the single one, and this is not surprising since now we must 
reconstruct four blocks of variables. By looking at the results 
reported in Table 6 we see that the two initializations produce 
in all cases two sequences of iterations converging to distinct 
results. Even if, in some cases, the two values of the objective 
function are very close, the corresponding points are definitely 
different, thus implying the existence of several minimizers or 
stationary points with very close values of the objective func¬ 
tions. 

It is interesting to remark that, while in the case of the bi¬ 
naries the best results are provided by the second initialization, 
now they are provided by the first one, based on the autocorre¬ 
lations of the ideal PSFs. The highest reconstruction errors are 
obtained in the case of the lowest SR, as one should expect. We 
also remark that in the case of the second initialization we used 
a larger number of iterations because the convergence is slower 
than in the case of the first initialization. From the comparison 


of the results obtained for the binaries with those obtained for 
the star clusters we deduce that the problem of the initial PSFs is 
essentially open; therefore, in the case of practical applications, 
one should try with different initializations, using also physical 
intuition in their choice. 

As a final comment, all the values of the objective function 
corresponding to the best solutions are higher than those corre¬ 
sponding to the other ones. 


5. Conclusions 


In this paper we extend to the case of Fizeau interferometry 
a blind deconvolution method previously proposed for single 
aperture telescopes and we validate the method in both cases, 
called respectively multiple image and single image case. 

It is well-known that the problem of blind deconvolution is 
extremely ill-posed and the introduction of constraints on PSF 
and object does not exclude the existence of several local min¬ 
ima, stationary points etc. In our approach the most signifi- 
ca nt constraint is the SR constr aint on the PSFs, as suggested 
in Desidera and Carbillel ( 20091) . This constraint excludes the 


trivial solution of a delta function for the PSF and image for the 
object. 

From our numerical analysis it turns out that the problem of 
Fizeau interferometry is more difficult than the problem of sin¬ 
gle aperture telescopes. The reason may be twofold. On one 
hand the number of variables to be reconstructed is larger and it 
is known that in the minimization of a nonconvex block-convex 
function the theoretical results are weaker when the number of 
blocks is greater than two. On the other hand the PSFs are very 
structured due to the presence of interference fringes so that 
if the initialization does not contain sufficient information on 
these structures it is difficult if not impossible to obtain accept¬ 
able results. 

An astonishing feature already observed in the single aper¬ 
ture case and confirmed in the present paper is that very often 
the value of the objective function corresponding to a sensible 
solution is greater than the value corresponding to an unaccept¬ 
able one. Obviously it is impossible to verify if these points 
correspond to local minima or to stationary points. In any case 
this result raises the issue if global minima, in case they were 
computable, provide sensible solutions or not. 

In summary the results of this paper open a large number of 
problems; however we think that the proposed method, which 
has a sound mathematical foundation, is very flexible and can 
help to investigate these problems. Moreover, last but not least, 
it can be extended to introduce regularization terms both on the 
object and on the PSF (or PSFs) thanks to the high flexibility 
of SGP, which is the basic tool in our approach. It is suffi- 
cient to mod i fy the scaling factor along the lines suggested in 
Lanteri et alJ (120021) . Obviously, in such a case, the additional 
problem arises of the choice of the regularization parameters. 

We conclude by remarking that the IDL routines implement¬ 
ing our method are available on request and will be included in 
the forthcoming versi on of the Software Package AIRY (v.6.1) 
( Carbillet et al. . 2014 ). 
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Table 6: Multiple image case - The reconstruction eiTors in the case of the two models of star cluster. The structure is similar to that of Table|5]but now we give the 
errors on the three PSFs and the normalized objective function is defined by 2Jo/3N^. 


Star Cluster 

SR 

Init. 

MARE 

RMSEqo 

RMSE60“ 

RMSEi2o» 

Jo norm. 

IT 


0.77 

A 

0.35% 

2.14% 

4.28% 

4.19% 

0.3049 

1000 

OC 

C 

0.59% 

4.16% 

7.67% 

7.62% 

0.2997 

5000 

0.46 

A 

0.81% 

3.43% 

6.07% 

6.00% 

0.1321 

2000 


C 

0.89% 

3.74% 

7.77% 

7.81% 

0.1237 

10000 


0.77 

A 

0.38% 

1.98% 

3.46% 

3.38% 

0.7597 

3000 

GC 

C 

0.71% 

11.00% 

11.06% 

12.97% 

0.7459 

5000 

0.46 

A 

1.04% 

5.63% 

9.69% 

9.38% 

0.3557 

6000 


C 

0.90% 

25.81% 

16.37% 

17.94% 

0.3043 

10000 
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